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Abstract 

The dynamics of an interacting Fermi gas of atoms at sufficiently high temperatures can be efficiently studied 
via a numerical simulation of the Boltzmann equation. In this work we describe in detail the setup we used 
'. recently to study the oscillations of two spin-polarised fermionic clouds in a trap. We focus here on the 
■ evaluation of interparticle interactions. We compare different ways of choosing the phase space coordinates 
of a pair of atoms after a successful collision and demonstrate that the exact microscopic setup has no 
■4-^ ■ influence on the macroscopic outcome. 

C ; 
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1. Introduction 

The Boltzmann simulation of the semiclassical regime of trapped Fermi gases [l[ has become a crucial 
tool to understand recent experiments 0, 0, 0, [1] where spin-polarised atomic clouds undergo strong col- 
lisions. The simulation is one of the very few methods which allow us to calculate accurately e.g. spin 
transport coefficients in ways that can then be quantitatively compared with experiment. In a previous 
work we studied the collision of two spin-polarised fermionic clouds and obtained excellent agreement 
with experimental results. The goal of the present paper is to describe in detail the numerical setup that 
was used there and to present tests of the simulation and other information which is potentially very useful 
for those wishing to implement this important method. 

Our setup is related to that of but we introduce several new modifications 'sj. For instance, we 
place an artificial grid over the continuous coordinate space that allows us to evaluate particle collisions 
efficiently. However, the main difference between this work and [7] is the method used to choose the new 
positions and momenta of two colliding particles after the collision. A quantum mechanical collision is a 
fundamentally random process with constraints given by the symmetries of the system. There are several 
ways to implement these constraints within the framework of a molecular dynamics simulation. We will 
discuss the physics of several collisional setups and also compare the numerical results associated with them 
in detail. We find that, although they exhibit very similar macroscopic behaviour, there are distinct technical 
advantages which make some of the methods preferable over others. 

In Sec. [5] of this paper we introduce the general numerical setup. We describe the method of test 
particles and explain what conditions need to be fulfilled for a semiclassical collision to take place and 
how to efficiently search the system for suitable collision partners. The different ways of choosing the new 
positions and momenta after a collision are discussed in detail in Sec. [31 In Sec. H] we show how to determine 
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the optimal simulation parameters and present several tests of the numerical method, in particular also of 
the different coUisional setups. A summary of our results and an outlook on future work can be found in 
Sec. [SI 



2. Numerical setup 

2.1. Definitions 

We consider a system of two-component fermions with equal mass m, labelled by the spin index s = {^. ^} 
and work in units in which ft = /cb = 1. Fermions of opposite spin can interact via s-wave scattering and 
the cross section is given by 

where a is the scattering length and Prci = Pt ~ Pi is the relative momentum of the two atoms. We assume 
that the system is in the normal phase and that the temperature is sufficiently high, so that the two spin 
distributions can be described semiclassically in terms of functions fs{r,p,t). Within the local density 
approximation the equilibrium distribution for fermions is given by the Fermi-Dirac distribution 

f(FD)/ X ^ I C2^ 

where T is the temperature and fig the chemical potential defined by the normalisation condition 

N. J d\n,{v,t) ^ Jd'rJ ^f,{r,p,t). (3) 

The atom number for each species is assumed to be equal, Nf = = N/2, and will be held fixed during 
the simulation. The trapping potential is assumed to be harmonic, 

Vir) - ^miulx' + uy + ulz^), (4) 

with the three trapping frequencies w^;, ujy and lUz. In the following we will only consider the isotropic case 
LOx — i-tJy — ijJz or a cigar shaped trap with < w^. = uiy = lj^. We denote the geometric average of the 
trap frequencies as loq = (w^w^w^)-'^/^. 

We define the Fermi energy from the energy distribution of the non-interacting gas at zero temperature, 
Ey ^ Ty ~ (3A'^)^/'^a;o. This energy scale marks the frontier between the classical and quantum regimes. 
For T > Tp the fermions behave essentially like classical particles and can be described by the Maxwell- 
Boltzmann distribution, 

/rB)(r,p) = Ar^||e-(pV2™+y(r))/T^ 

The Fermi energy can also be used to determine the typical scales of the cloud, ep — fc|/2m = mujfRf/2, 
where fcp is the Fermi momentum and Ri ~ kY/{muJi) are the Thomas- Fermi radii in the three spatial 
directions. These quantities give the widths of the zero-temperature momentum and density distributions 
respectively and are therefore useful to describe the extent of the cloud in momentum and coordinate space. 
The time-evolution of the distribution function /s(r, p, t) is given by the Boltzmann equation, 

dtfs + (p/m) • Vrfs - VrV ■ Vpfs - -nfsjsi (6) 

where the left-hand side represents the propagation of the atoms in the potential and the right-hand side 
stands for the collision integral, which depends on the particle statistics. Here, we consider only collisions 
between atoms carrying opposite spins. Indeed, for ultracold fermions, the Pauli principle forbids interactions 
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between particles of the same spin fT]. At temperatures T >Tf when the Maxwell-Boltzmann distribution 
is applicable, the collision integral takes the form 

Iclassifs, fs] = J J TO^"^ [fsfs- f'sfir], (7) 

where the indices s and s label the two colliding atoms of opposite spin, the primed variables refer to 
quantities after the collision and is the angle between the incoming and outgoing relative momenta. The 
Pauli exclusion principle only allows fermions to scatter into a previously unoccupied quantum state. This 
reduces the scattering probability by the Pauli blocking term proportional to (1 — /^)(1 — /^). Taking this 
into account, the collision integral reads 

At temperatures T > Tp the Pauli terms are close to one, and therefore the collision integral (|8]) tends to 
the classical form ([7]). 

2.2. Test particles 

To efficiently simulate the evolution of the continuous distribution function /^(r, p, t) we use the method 
of test particles. These are point-like particles which form a discrete approximation to fs{r,p,t) through 
(5-functions. In order for this approximation to be accurate we will represent each fermion by several test 
particles 0, The higher the ratio N/N of test particles to atoms, the more precisely the continuous 

distribution function will be approximated, 

N 

/,(r, p, t) ^ ^(2^)35(r - r,{t))S{p p. (i)). (9) 

Physical observables need to be rescaled, for instance the test particle cross section becomes a = a{N/N). 
A generic thermal expectation value of a single-particle observable X(r, p) can be easily calculated within 
the test particle picture, as the integration reduces to a sum over all test particles, 

W = H^J d'r^,fs{r,V,t)X{v^v) (10) 



1 ^ 

^^X(r„p,). (11) 



i=i 

We introduce a discrete time step Ai, such that during each time step the test particles propagate 
without colliding, following their classical trajectories. At the end of each time step collisions between them 
are evaluated. In a harmonic potential the trajectories are given by 

ri(tn+i) = ri{tn)cos{LLiiAt) + {pi{t„)/muji)sin{ujiAt), (12) 
Pi{tn+i) = Pi{tn) cos{ujiAt) - ri{tn)mujiSm{ujiAt). (13) 

Note that since the time step is fixed, the trigonometric functions only need to be evaluated once during 
the entire simulation, so that using the exact solution is more efficient than using the Verlet algorithm, as 
for instance in Refs. 7, 9], in which case the accelerations need to be recalculated for every time step. The 
Verlet algorithm is more general as it is applicable for any potential. But since in this work we will only 
consider harmonic potentials, we will use the exact solution instead. 

We evaluate the probability of a two-particle collision in the same way as described in fzj . First we must 
test whether a given pair of test particles reaches the point of closest approach during the present time 
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step. This condition is important to prevent particles from attempting to collide with each other repeatedly 
over several consecutive time steps, an issue which will be further addressed below. If the closest approach 
condition is true we check if the minimal distance dmin of the test particles fulfils the classical condition for 
scattering: Trd^j^ < a. If this condition is also satisfied we propose a collision at the time of closest approach. 
However due to Pauli statistics, even if the classical conditions for scattering are fulfilled, a collision can 
only take place if the new state of the particles was previously unoccupied. To take this into account we 
calculate the quantum mechanical scattering probability given by the Pauli term (1 — — /-^) and accept 
or reject the collision according to this probability. Clearly, the point-like particle picture is unsuitable for 
the calculation of this probability. To return to a continuous distribution we therefore have to smear out 
the (5-functions representing the test particles, e.g. by Gaussians in position and momentum space: 

e-(p-P.)V»p 

S{p-P^) ^ . ^ .3 (14a) 

f,-{^-x,f/wl „-{y-yi?lwl -(z-z,f/wl 
6{r-r,) -> y= y= . (14b) 

^/■7^Wx \jT^Wy ^/-KWz 

The widths of these Gaussians, Wx, Wy, Wz and Wp, need to be tuned so that on the one hand fluctuations 
due to the discrete nature of the test particle picture are smoothed out, but on the other hand the physical 
structure of the distribution function fs remains preserved 01 • The first condition is equivalent to 

WpWr > {N/Ny/^, (15) 

where we introduced Wr — {wxWy'WzY^'^ , the geometric average of the spatial widths. The second condition 
implies Wi ^ Ri and Wp ^ kp. 

We also require that the smearing by the Gaussian functions preserve temperature-dependent degeneracy 
effects of the Fermi distribution, particularly close to the Fermi surface. This implies 

Wp ^kFiT/TF) and Wi <s: RiiT/fp). (16) 

Note that the smearing width in momentum space is isotropic, while in position space the smearing width can 
be different depending on the spatial direction, if the corresponding trap frequencies are unequal. Since the 
Thomas-Fermi radii are inversely proportional to the corresponding trap frequencies it is sensible to choose 
Wi — WrLdo/uJi for the spatial widths. Furthermore Eq. (IT5|) together with the definition of the Thomas- 
Fermi radii imply Wp = mLOoWr- Hence all four smearing widths can be reduced to only one free parameter. 
At very low temperatures the margin given by the conditions (1151) and (jl6p becomes so narrow that it is 
impossible to find smearing widths satisfying both conditions, without having to significantly increase the 
number of test particles. This limits the applicability of this setup to temperatures above approximately 
0.2Tf for N/N = 10. Moreover, at very low temperatures the system undergoes a phase transition into a 
superfiuid state. This method does not include the relevant degrees of freedom of the superfluid Fermi gas 
and is only applicable in the normal phase. 



2.3. Auxiliary grid 

The main numerical challenge is to efficiently evaluate collisions between the test particles. The total 
number of pairs of opposite spin is 7V^/4, which is an unfavourable scaling given that we want to use large 
particle numbers and a high test particle to particle ratio. In this work we develop a more efficient method 
than to check all possible test particle pairs. The key observation is that since the cross section is decreasing 
with increasing relative momentum it can never exceed (imax = Ana^N/N and consequently the maximal 

distance two colliding test particles can have is dmax — 2ay^ N/N. Having this in mind we superpose a cubic 
grid with cell size d^ax on the continuous space. The grid has finite extent which can be set by demanding 
that a certain proportion, for instance at least 95%, of the test particles are within the grid. For a cigar 
shaped trap the grid must have larger extent in the axial direction. At the end of each time step we move 
systematically through all grid cells starting in one of the corners and note all possible collision partners 
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that fulfil the classical scattering conditions. To make sure that we do not miss collisions due to boundary 
effects, for each particle we check not only all opposite spin particles in the same grid cell, but also in all 
neighbouring cells (the ones sharing a face, an edge or a vertex with the given cell). This ensures that all 
particles in a sphere of radius dmax around a given particle are definitely accounted for. This makes a total 
of 3^ — 27 cells for each particle, however to avoid double counting we only need to evaluate cells in the 
positive direction, which means on average 14 cells per particle. 

A small systematic error source remains with this setup. If the relative velocity of two particles is large 
they can be in non-neighbouring cells at the beginning and at the end of a time step, although in the course 
of the time step they come within each other's allowed collision range. Such a possible collision will then 
not be accounted for. However this systematic error can be minimised by choosing the time step to be 
sufficiently small and also by choosing the cell size to be larger than dmax- Also note that for large relative 
velocities the cross section is small and collisions between very fast particles are rare events. 

After having searched the entire grid for classically allowed collision pairs we proceed to choose which 
collisions will indeed take place. To do so we consecutively select random pairs from the list. We then prop- 
agate both particles to the point of their closest approach, let them scatter (the exact setup for determining 
the new positions and momenta after scattering is described in Sec. [3]) and then propagate them back to the 
original time. To account for quantum statistics we then calculate the Pauli blocking factors using the new 
distributions = /s(r',p'). The continuous distributions are obtained by replacing the (5-functions in ^ 
according to (fH)) . If we obtain a value > 1 from the summation we set f!, = l. The probability that the 
collision is accepted is then given by {1 — f!,){l — f-^). Regardless if a collision is accepted or rejected neither 
of the particles concerned is allowed to collide again with another particle during the present time step. If 
the collision is accepted we keep the new positions and momenta. If the collision is rejected we return to 
the values before the collision. This procedure is repeated until all possible pairs have been evaluated. 

3. Collisional setup 

Since the picture of colliding point-like particles with well-defined positions and momenta is a classical 
interpretation of a quantum mechanical scattering process, it is unsurprising that there must be some 
ambiguity in the implementation of the collisional setup. In the semiclassical particle picture each collision 
has 12 degrees of freedom: three position and three momentum components for each of the two particles, or 
equivalently three components of the total and relative positions and momenta. In quantum mechanics we 
consider wave packets rather than particles and concepts like particle position or momentum are not well- 
defined. Instead the system is described by operators which, depending on the symmetries of the system, 
commute with the Hamiltonian and define quantum numbers corresponding to conserved quantities. The 
concept of a trajectory does not exist, as a particle is not localised in phase space but rather smeared out 
over a certain phase space volume in accordance with Heisenberg's uncertainty principle. Therefore, in the 
collisional setup, it is not necessary to preserve for instance the exact positions of the two atoms during 
a collision. In the presence of an axially symmetric external potential for instance, the only conserved 
quantities are the total energy and the angular momentum component L^. 

3.1. Angular momentum preserving setup 

The method used in this work and in \6\ can be motivated as follows. As collisions are local, we can 
disregard for the moment the external potential and go to the centre of mass frame of the two particles. 
Motivated by the analogue of classical scattering we wish to conserve the total momentum, the total energy 
and the total angular momentum of the system during the collision. Conservation of the total momentum 
p = Ps 4- Ps and the total energy E = E^kin — {pi +p'^)/2m together imply the conservation of the modulus 
of the relative momentum p^c\ = |Ps — Ps|, since p^ +Prc\ — 2(Ps +-P|) — 4:mE. The direction of the relative 
momentum vector can change during the collision. Finally we must also conserve the angular momentum 
L = Trei X Prci- Wc Can Satisfy all these constraints by rotating both the relative momentum and relative 
position vectors by the same arbitrary angle around the L-axis. The angle of this rotation is the only degree 
of freedom of the collision and is determined at random. From the new values pj.^j and rj.^j the new values of 
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the momenta and positions of the individual particles can be recovered using total momentum conservation 
and the centre of mass coordinates respectively. 

So far we have ignored the external potential, which would be justified if the collisions took place when 
the atoms were infinitesimally close to each other and therefore experiencing the same potential. However, 
the two colliding particles are not exactly at the same position before the collision and their relative position 
changes after a successful collision. Thus potential energy is not conserved and total energy conservation 
is not exact during a collision in general. Nevertheless, we will show below that such changes in the total 
energy are negligibly small for almost any collision. Furthermore, these small changes cancel each other in 
a many-particle system which experiences many collisions. 

3.2. Energy preserving setup 

It is possible to preserve energy conservation exactly by employing a different setup, for instance as in 
. In this reference, the relative position stays fixed during a collision and the relative momentum vector is 
rotated by a random angle in space (such a rotation has two degrees of freedom). As a direct consequence 
of the unrestricted rotation this setup violates angular momentum conservation. While it is true that this 
violation would occur naturally in non-axially symmetric potentials even for non-interacting atoms, this 
coUisional setup provides an additional, unphysical change of angular momentum. At any rate, since in this 
work we consider either isotropic or axially symmetric potentials, either the total angular momentum or its 
axial component are conserved. 

4. Tests and optimal parameters 

How accurately the numerical setup represents the physical picture depends crucially on the values of 
the simulation parameters, in particular N/N, the time step At and the smearing widths Wr and Wp. For 
all tests described below we use N = 10000 atoms and N/N — 10, which is sufficient for temperatures 
T > 0.2Tf. 

The optimal value of At depends on the physical parameters of the system. Obvious requirements are 
that the time step must be smaller than the typical time between two collisions and that the average distance 
travelled during a time step must be much smaller than the diameter of the cross section, (uioi) At <C y^Jaj/n. 
Another constraint is that the time step must be smaller than the half-period with respect to the largest 
trap frequency, At < Tr/wmax, but unless the aspect ratio is extremely large this condition is much weaker 
than the other ones. There is no lower bound on the time step, however the simulation slows down with 
decreasing At. 

We first describe the tuning of the parameters and the corresponding tests of the simulation with the 
angular momentum preserving collisional setup 13.11 We perform the same tests as in [7] to demonstrate 
that even with the different collisional setup we obtain very good agreement. Then we explicitly compare 
the different collisional setups in Sec. 14.41 

4.1. Equilibrium collision rates 

To obtain the correct dynamical properties, for instance the damping time of excitation modes, we need 
to ensure that the equilibrium collision rate observed in the simulation corresponds to the correct theoretical 
value. The number of collisions in a given time interval can be very easily obtained from the simulation, 
simply by counting all test particle collisions and then multiplying by the ratio (N/N). The theoretical 
value for the collision rate 7 in the presence of Pauli blocking is given by the following integral. 

After inserting the Fermi-Dirac distribution this integral can be calculated numerically 7]. The numerical 
setup also provides a powerful tool to artificially switch off Pauli blocking. This allows us to separately 
check for errors related to the general numerical setup and the calculation of the Pauli blocking factors. 
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Without Pauli blocking the integral is simpler and can be solved analytically for the Maxwell-Boltzniann 
distribution, 

= N,N,^(l + ^e^EiL^)), (19) 

where Ei(a;) = J^^{e*' /t)dt is the exponential integral. Furthermore we can obtain the theoretical prediction 
for the Pauli blocking probability by solving the integral ([T5| for the Fermi-Dirac distribution. The proba- 
bility ppauii that two classically colliding fermions will indeed scatter is then given by the ratio of 7biock to 
this integral. 

To find the optimal value for At for each system we measure the collision rate in the absence of Pauli 
blocking for decreasing values of the time step and compare it to the theoretical prediction. The time step is 
small enough when we reach good agreement. It is important to switch off Pauli blocking for the tuning of 
the time step, as in the presence of Pauli blocking the collision rate is sensitive to the values of the smearing 
widths, which can obscure inaccuracies due to a time step that is too large. After having found the optimal 
time step we check that repeated (unphysical) collisions of the same particle pair are rare. This has always 
been found to be the case with our coUisional setup. We then use this optimal value for At to establish the 
optimal values of the smearing widths. We first identify the allowed interval for Wp and Wr given by the 
conditions ([T5|) and ([T5| and perform the same collision rate matching as described above, this time with 
Pauli blocking switched on. We find that the optimal widths lie between Wr = Wp / {mLOo) = l.O^ho for the 
lowest and Wr = Wp/lmujo) — 2.0/ho for the highest temperatures used in our analysis, where Zho — 1 / ^mioii 
is the natural harmonic oscillator length unit. 

The measured collision rates for the optimal choice of parameters with and without Pauli blocking to- 
gether with the theoretical predictions are shown in Fig. [T] For sufficiently high temperatures the agreement 
is very good. At very low temperatures it becomes increasingly difficult to resolve the conditions (|T5|) and 
(1161) for the Gaussian smearing widths simultaneously. For larger values of the scatteringlength this problem 
gets worse since test particles with a large relative distance can scatter on each other [7| and therefore the 
continuous distribution function needs to be resolved accurately over larger scales. 




Figure 1: (Colour online) The equilibrium collision rates per particle (in units of (Jq) with and without Pauli blocking, as well 
as the Pauli probability for a successful scattering versus temperature for Ifepa] = 1 (loft) and |fcpa| = 2 (right). The lines 
correspond to the theoretical prediction and the symbols to the values obtained with the simulation. 



^.2. Equilibrium energy distributions 

Another important test is to check that the system thermalises to the correct equilibrium energy distri- 
bution, independently of the initial distribution. In the presence of Pauli blocking the energy is distributed 
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according to the Fermi-Dirac distribution, whereas without Pauh blocking the particles will be distributed 
according to the Maxwell-Boltzmann distribution. Figures [5] and |3] show the results of this test for a low 
temperature system with \kpa\ = 1 and isotropic trap frequencies. The parameter values for the time step 
and the smearing widths are optimal. We performed two tests of the thermalisation. First we initialised 
the system according to the Fermi-Dirac distribution for T — 0.2Tf (Fig. and ran the simulation without 
Pauli blocking. After a short time the system thermalised according to the Maxwell-Boltzmann distribution 
for T = O.SlTp. Note that the temperatures of the two distributions are not necessarily equal, since the 
equipartition theorem does not hold for the Fermi-Dirac distribution. When changing from one distribution 
to the other the average energy of the system remains conserved and hence the temperature of the new 
equilibrium state is different. Figure [3] shows the corresponding results for the reverse situation: the initial 
distribution is Maxwell-Boltzmann with T = O.SITf and Pauli blocking is switched on. It is clearly visible 
from both figures that the correct equilibrium distribution is always attained at the end of the simulation. 
This agreement improves further at higher temperatures. 





Figure 2: (Colour online) The equilibrium energy distributions without Pauli blocking. The top panel shows the energy 
distributions on a linear scale, the bottom panel shows the energy distributions scaled by lJq/E'^ on a logarithmic scale. The 
start distribution (left) is Fermi-Dirac (F.-D.) at T = 0.2Tp and the end distribution (right) is Maxwell-Boltzmann (M.-B.) 
with the same average energy but at T = 0.31Tp, as expected. 



4-3. Collective excitations 

Collective excitations emerge when a many-particle system is perturbed away from equilibrium. Here 
we will discuss three different excitation modes, the sloshing mode (also known as dipole or Kohn mode), 
the breathing mode (monopole mode) and the quadrupole mode. We will confirm that the simulation gives 



8 





Figure 3: (Colour online) The equilibrium energy distributions with Pauli blocking. The top panel shows the energy distributions 
on a linear scale, the bottom panel shows the energy distributions scaled by uj^/E'^ on a logarithmic scale. The start distribution 
(left) is Maxwell-Boltzmann (M.-B.) at T = O.SlTp and the end distribution (right) is Fcrmi-Dirac (F.-D.) with the same average 
energy but at T = 0.2Tp, as expected. 
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Figure 4: Simulation of the equilibrium sloshing mode {z)/lz, where Iz = 1 / ^mu)z ■ The mode is undamped and the oscillation 
frequency equals Loa- 
the correct frequencies and damping properties of these modes. The fohowing tests were performed for a 
spherical trap ujx = tOy = uJz = <-^o- 

The sloshing mode is excited by a small displacement of the centre of mass from its equilibrium position, 
or equivalently by a short-lived force represented by an additional linear term in the potential. The time 
evolution of each of the three centre of mass coordinates (r^) is well-known to be an undamped oscillation 
with the corresponding harmonic oscillator frequency . Figure |3] shows such an oscillation for a system at 
Ifcpal = 1 and T 0.2fF. 

The breathing mode can be excited by compressing or expanding the cloud. In a spherical trap this yields 
an undamped oscillation of (r^) with frequency 2aJo around its equilibrium value Figure [5] illustrates 
this mode for a system with jfepil — 1 and T = 0.2Tf. 

By perturbing the system via a short-lived small increase in one or several of the trap frequencies we 
can excite the quadrupole mode. We excite the quadrupole mode Q{t) — (x^) — (j/^) by applying the 
perturbation l^px = —cx and Apj, = cy with c = 0.2ma;o in the same way as in ;7:]. Unlike the sloshing 
and the breathing mode this mode is damped. The frequency of the quadrupole mode at high temperatures 
approaches the ideal gas value 2ajo, while at low temperatures it is closer to the hydrodynamic frequency 
a/2wo- Figure [5] shows the quadrupole mode for \k-pa\ = 1 in the high and in the low temperature regime. 
In both cases the damping of the mode is clearly visible. The frequency of the oscillation can be extracted 
from the corresponding Fourier transform (5(w) = (iiQ(i)e'"* and is in agreement with the theoretical 
prediction while the damping can be extracted from its imaginary part. We see that the damping is stronger 
at the lower temperature T = 0.4Tf when collisions are more frequent (see Fig. [l]for a plot of the collision 
rate). 



^.4-. Comparison of different collisional setups 

To study the impact of the collisional setup on the outcome of the simulation we have implemented the 
setups 13.11 (angular momentum preserving) and 13.21 (energy preserving) and confirmed that they generate 
the same results for both equilibrium and non-equilibrium systems. First we analyse the magnitude of the 
changes in the total energy with the angular momentum conserving setup. Figure [7] shows a typical histogram 
of the relative change of the energy of a particle pair AE/Einit = (^^finai — Einit)/Einit after a successful 
collision. The histogram is sharply peaked around zero. In this typical example the majority of collisions 
(approximately 84%) conserve the energy of the colliding pair with an accuracy of up to \AE/Einit \ < 10"''. 
Since collisions are frequent and the energy changes have random sign we also observe large cancellation 
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Figure 5: Simulation of the normalised equilibrium breathing mode (r^)/(r^)eq. The breathing mode is undamped and the 
oscillation frequency equals 2ujz ■ 




Figure 6: Left: Simulation of the normalised equilibrium quadrupole mode (3(t)/(r^}eq at high temperature T = IT-p (top 
panel) and at low temperature T = 0.4Tp (bottom panel). Right: The imaginary part of the corresponding Fourier transforms 
of the numerical data gives information about the oscillation frequency and the damping. 
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effects, such that the total energy of the system is conserved with an accuracy of the order of 10 ^ at any 
given time. This is smaller than for instance the energy deviations due to Verlet's algorithm quoted in [?[. 

40000 
35000 
30000 
25000 
20000 
15000 
10000 
5000 


-0.0001 -0.00005 0.00005 0.0001 

AE/Ei,i, 

Figure 7: Probability density (obtained from the corresponding normalised histogram) for the relative change of the energy of 
a colliding particle pair AE/Sinit = (Bflnal - £^init)/Einit for Ifcpal = 1, T = 0.4Tp, wx/i^z = 100 and TV = lOA^ = 100000. 

To understand the good accuracy of energy conservation we also make a qualitative theoretical estimate. 
During the collision the kinetic energy is conserved, such that it is sufficient to consider the potential energy 
V of the two particles, 

^ - y [^1(4 + yl + 4+ vD + + ^D] , (20) 

where {xi,yi,Zi) are the coordinates of the particle i and lu± = uj^ = ojy is the radial trap frequency. We 
can recast this expression in terms of the centre of mass and the relative coordinates {X, Y, Z) and {x, y, z) 
respectively, 

V = m [u^liX^ + Y') + ^IZ^] + ^ + ,/) + ^Iz^] . (21) 

The position of the centre of mass is not affected by the collision. The variation of energy is thus only due 
to the contribution of the relative motion Vrci defined by 

y,ei = jK(x2+Z/^)+^^^] (22) 

We now define a second frame (a;', y', z') such that the z'-axis is parallel to the relative angular momentum L 
of the pair. In addition, since the trap is rotationally invariant around z, we can choose the initial coordinate 
system such that L is in the (x, z) plane, in which case y' — y. If ^ is the angle between L and the z axis, 
the two sets of coordinates are related by the following relations 

X = x' cos 6 + z' sin 9 

y ^ y' 

z ~ z cos9~x SU10. 

By conservation of angular momentum, the relative coordinates evolve in the {x' ,y') plane. We thus have 
z' = 0, cc' = dcosip and y' = d sin 1^9, where d is the relative distance between the two particles. In terms of 
d and ip, the relative potential energy is 

Kci = ^ |-^-(u;f-wi)sin^6'cos2v'. (23) 
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During the collision d and 9 are conserved. The variation of energy is thus only due to the rotation in the 
{x',y') and we obtain 

AE = ^^(w^ - wi) sin^ 9 [cos^ (pf - cos^ ipi] , (24) 

where tpi and iff correspond to the initial and final values of ip respectively. As an approximation, we 
normalise to the average total energy of the pair, E w 6T, and obtain 



As shown above the distance between the particles can be at most dmax = '^o.y N/N. Hence for the range of 
parameters considered in this study the constant prefactor is of order < 10~^. The trigonometric functions 
further reduce the variation of the energy and are responsible for the pointed shape of the histogram. Note 
that the method of test particles is important for a small energy variation. Due to N being larger than N 
the cross section and hence the average distance between two scattering partners is reduced, which leads 
to a smaller AE/Einit (in our case by a factor of N/N = 10). To demonstrate this effect Fig. [8] shows the 
variance of the energy change as a function of the ratio N/N. The data is well-fitted by a line with slope 
—2 on a double-logarithmic scale, which confirms the relation ([25)) . Analogous scaling laws can be shown 
also for the dependence on the temperature and the scattering length. 

10-7 , . . . ^ 



in 

a 
> 



10 

N/N 

Figure 8: Double-logarithmic plot of the variance of the energy change during a collision Var(Ai?/i?init) for |fcpa| = 1, T = Tp 
and ujj_/lJz = 100 as a function of the test particle to particle ratio N/N. The line corresponds to a linear fit with constant 
slope —2 on the double-logarithmic scale. 

We also compared the time evolutions of the centre of mass coordinates in response to a perturbation. 
The differences between the two setups were found to be smaller than the statistical fluctuations. Taking 
these results into account we conclude that the two setups 13.11 and 13.21 are equivalent on a macroscopical 
scale and we have the freedom to choose whichever option appears more convenient. 

We adopt the angular momentum conserving setup 13.11 not only because of this property, but also since 
the setup from ^ has a small technical disadvantage. Since the new direction of the relative momentum 
vector is chosen uniformly on a sphere, in many cases the particles are found to approach each other again 
after a successful collision. This implies that in the following time step they are likely to undergo another 
collision. The proportion of such events compared to the total number of collisions can reach up to 40% 
at the lowest temperature considered in this study and amounts to around 10% at T = Tp. To avoid this 
overcounting of collisions one needs to implement an additional routine that prevents particles from colliding 
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with each other repeatedly within short time intervals. In other words the molecular chaos assumption does 
not hold with the setup from and needs to be enforced artificially. In our setup repeated scattering is 
so rare that this effect can be neglected. This is because the relative position and the relative momentum 
vectors of the particle pair are rotated by the same random angle and hence their relative angle is preserved. 
After a collision takes place the particles continue to move on trajectories that diverge from each other. 

4-5. Outlook 

In addition to the two collisional setups described above we propose a third very general coUisional setup 
which can be explored in future work. 

To preserve macroscopical averages we keep the centre of mass coordinates and momentum fixed and 
work with the six degrees of freedom for the relative momentum and position. Based on the idea of a 
delocalised particle pair, we assume that the relative phase space coordinates of the particles are distributed 
according to a Gaussian. We then choose the new relative position and momentum randomly according 
to the Gaussian probability distribution with the additional constraints given by the symmetries of the 
problem. Energy and Lz conservation for instance reduce the problem by two degrees of freedom. If energy 
and total angular momentum are conserved the problem is reduced by four degrees of freedom and so on. 

It is important to stress however, that the task of defining a probability measure on a complicated 
hypersurface in the phase space (in our case defined by the constraints of energy and angular momentum 
conservation) is a mathematically highly nontrivial problem and will pose a challenge for the numerical 
implementation. 

5. Conclusion 

We have implemented a Boltzmann equation approach to the simulation of the semiclassical regime of 
two-component Fermi gases. We described how to implement the method using test particles and stressed 
the advantages of the auxiliary grid. We also presented an extensive range of tests for the simulation. The 
main focus of this work was a discussion of the different possibilities for implementing collisions, namely 
the angular momentum conserving and the energy conserving setups. We studied in detail the error in 
the energy of the former scheme and showed that it is indeed extremely small, so that it becomes a valid 
alternative way of calculating the collisional integral and has distinct advantages since it minimises repeated 
scattering events. In future work we intend to apply the method to systems with components of unequal 
mass which bring a greater variety of dynamical regimes and which can provide a more stringent test of the 
Boltzmann equation approach to these systems. 
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